#ifndef SPECIES_H
#define SPECIES_H

#include <vector>
#include <string>

#include "CellParticles.h"
#include "PartBoundCond.h"
#include "../Input/PicParams.h"
#include "../Input/PicParams.h"
#include "../Parallel/ParallelMPI.h"
#include "../PhysicalField/PhysicalField.h"
#include "../Field/Field.h"
#include "../Field/Field3D.h"
#include "../Tool/Log.h"

class PhysicalField;
class Pusher;
class Interpolator;
class Projector;


class Species
{
public:
    //Species creator
    Species(PicParams* params_, ParallelMPI* pmpi_, int ispec);

    //Species destructor
    virtual ~Species();

    //calculate the particle dynamics (interpolation, pusher, projection)
    //only caculate the number density, no electric currents
    virtual void move(PhysicalField* fields){};

    //decide if particles hit boundary, cross MPI boundary, to move to another cell, or still stay in the cell
    virtual void sort_particles(){};

    //project particles to grid points, to calculate charge density
    virtual void project(PhysicalField* fields){};

    //After MPI exchange particles, delete useless particles
    virtual void correct_particles(){};

    virtual int get_number_of_particles()
    {
        return 0;
    }

    void get_maxwell_velocity(double velocity[3])
    {
        double vt = sqrt(2.0 * species_param.temperature * params->const_e / species_param.mass);
        double r1;
        double r2;
        for(int i = 0; i < 3; i++)
        {
            do{ r1 = (double)rand() / RAND_MAX; } while(r1 == 0.0);

            r2 = (double)rand() / RAND_MAX;
            velocity[i] = vt * sqrt( -log(r1) ) * sin(2.0 * M_PI * r2) + species_param.mean_velocity[i];
        }
    }

    PicParams* params;

    ParallelMPI* pmpi;

    int species_number;

    SpeciesParameter species_param;



private:


};

#endif
